图像的噪声抑制

参考文档

所谓图像噪声，是图像在摄取时或是传输时所受到的随机干扰信号。常见的有椒盐噪声和高斯噪声

椒盐噪声的特征：出现位置是随机的，但噪声的幅值是基本相同的^；高斯噪声的特征：出现位置是一定的（每一个点上），但噪声的幅值是随机的，且服从高斯分布^

图像噪声的抑制办法是设计噪声抑制滤波器，在尽可能保持原图信息的基础上，抑制噪声

  • 均值滤波器
  • 中值滤波器
  • 边界保持类滤波器

均值滤波器

对图像上的某一待处理像素给定一个“模板”，该模板包括了其周围的邻近像素，将模板中的全体像素的均值来替代此待处理像素 对图像上的每一个像素点执行此操作，这很容易通过一个卷积完成，其中卷积核为$\frac{1}{9}\begin{bmatrix}1 & 1 & 1 \\ 1& 1 & 1 \\ 1 & 1 & 1 \end{bmatrix}$

卷积操作

2D卷积 多通道版本

卷积核（滤波器）的每一层在各自的输入通道上「滑动」，产生各自的计算结果。某些通道的内核可能比其他通道的内核具有更大的权重，以便比强调某些输入通道（例如，滤波器的红色通道卷积核可能比其他通道的卷积核有更大的权重，因此，对红色通道特征的反应要强于其他通道） 然后将每个通道处理的结果汇在一起形成一个通道，整个滤波器产生一个总的输出通道 最后一个术语：偏置。偏置在这里的作用是对每个输出滤波器增加偏置项以便产生最终输出通道 卷积运算实现见./DB/data/conv.py

In [1]:
import cv2
import matplotlib.pyplot as plt
import numpy as np
In [2]:
#一个带噪声的图像示例
lena_noise=cv2.imread('./DB/image/11.jpg')
plt.figure(figsize=(8,8))
plt.imshow(lena_noise[...,::-1])
Out[2]:
<Figure size 576x576 with 0 Axes>
Out[2]:
<matplotlib.image.AxesImage at 0x19a58a99160>
In [3]:
import sys
sys.path.append('./DB/data/')
from conv import conv,media_filter,knn_smooth,min_gray_var,sigma_filter,gauss_kernel
In [6]:
t=conv(np.rollaxis(lena_noise,2),np.ones((3,3)),ratio=None)
plt.figure(figsize=(8,8))
plt.imshow(t,cmap='gray')
Out[6]:
<Figure size 576x576 with 0 Axes>
Out[6]:
<matplotlib.image.AxesImage at 0x22db5e90a20>

加权的均值滤波

均值滤波的缺点是，会使图像变得模糊，原因是它对所有点同等对待，在将噪声点分摊的同时，将景物的边界点也分摊了，为了改善效果，可以采用加权平均的方式来构造滤波器。几个典型的加权平均滤波器： image.png 之所以要把整数提出来，除了方便展示，也有计算上的原因，理论上整数模板比浮点数模板的计算更快

In [25]:
h=np.array([[1,2,1],[2,4,2],[1,2,1]])
t=conv(np.rollaxis(lena_noise,2),h,ratio=1/16)
plt.figure(figsize=(8,8))
plt.imshow(t,cmap='gray')
Out[25]:
<Figure size 576x576 with 0 Axes>
Out[25]:
<matplotlib.image.AxesImage at 0x1fc2f633d68>

好吧，我并未看到任何提升，尼玛感觉还下降了。。

中值滤波器

改进的加权均值滤波器效果提升并不明显，因此必须改换思路，中值滤波器就是一种比较有效的办法

因为噪声（如椒盐噪声）的出现，使该点像素比周围的像素亮（暗）很多。如果在某个“模板”中，对像素进行由大到小排列，那么最亮或最暗的点一定被排在两端，取排在中间位置上的像素点替换待处理像素的值，就可以达到滤除噪声的目的 image.png 对于图像，仍使用3*3尺寸的滑动窗口，步长为1

In [16]:
t=media_filter(np.rollaxis(lena_noise,2))
plt.figure(figsize=(8,8))
plt.imshow(np.rollaxis(t,0,3)[...,::-1])
Out[16]:
<Figure size 576x576 with 0 Axes>
Out[16]:
<matplotlib.image.AxesImage at 0x254b0a73a20>

效果确实不错，我试着多执行几次中值滤波，猜测结果会更好

In [18]:
t=media_filter(media_filter(media_filter(np.rollaxis(lena_noise,2)))) #执行了三次中值滤波
plt.figure(figsize=(8,8))
plt.imshow(np.rollaxis(t,0,3)[...,::-1])
Out[18]:
<Figure size 576x576 with 0 Axes>
Out[18]:
<matplotlib.image.AxesImage at 0x254aeac0630>
In [31]:
#再看一个椒盐噪声更严重的示例
t=cv2.imread('./DB/image/19.jpg')
plt.figure(figsize=(10,10))
plt.subplot(121)
plt.imshow(t[...,::-1])
t=media_filter(media_filter(media_filter(media_filter(np.rollaxis(t,2)))))
plt.subplot(122)
plt.imshow(np.rollaxis(t,0,3)[...,::-1])
Out[31]:
<Figure size 720x720 with 0 Axes>
Out[31]:
<matplotlib.axes._subplots.AxesSubplot at 0x254b3ed5be0>
Out[31]:
<matplotlib.image.AxesImage at 0x254b36fada0>
Out[31]:
<matplotlib.axes._subplots.AxesSubplot at 0x254b3d74128>
Out[31]:
<matplotlib.image.AxesImage at 0x254b3d7fe48>

注1：对于椒盐噪声，中值滤波效果比均值滤波效果好，两者比较： image.png

注2：对于高斯噪声，均值滤波效果比中值滤波效果好，两者比较： image.png

In [48]:
#实验：噪声图像来源见见参考1
t=cv2.imread('./DB/image/20.jpg')
plt.figure(figsize=(10,10))
plt.subplot(121)
x=conv(np.rollaxis(t,2),np.ones((3,3)),ratio=None) #均值滤波
plt.imshow(x,cmap='gray')
t=media_filter(np.rollaxis(t,2)) #中值滤波
plt.subplot(122)
plt.imshow(np.rollaxis(t,0,3)[...,::-1])
Out[48]:
<Figure size 720x720 with 0 Axes>
Out[48]:
<matplotlib.axes._subplots.AxesSubplot at 0x254be1caa20>
Out[48]:
<matplotlib.image.AxesImage at 0x254be1ff470>
Out[48]:
<matplotlib.axes._subplots.AxesSubplot at 0x254be1ff4a8>
Out[48]:
<matplotlib.image.AxesImage at 0x254be1ff3c8>

好吧，我完全没有看出对于高斯噪声均值滤波比中值滤波好哪怕一丁点？

边界保持类平滑滤波器

经过平滑滤波处理后，图像会变得模糊，分析原因，在图像上的景物之所以可以辨认清楚是因为目标物之间存在边界，而边界点与噪声点有一个共同的特点是，都具有灰度的跃变特性，所以平滑处理会同时将边界也处理掉

为了解决图像模糊问题，一个自然的想法是，在进行平滑处理时，首先判别当前像素是否为边界上的点，如果是，则不进行平滑处理，如果不是，则进行平滑处理

K近邻平滑滤波器

image.png 在模板中，分别选出5个与点1和点2灰度值最接近的点进行均值计算，这样就不会出现两个区域信息的混叠平均，因而可以达到边界保持的目的。具体的操作就是，以待处理像素为中心，做一个m*m的作用模板（m的值应为奇数），在模板中，选择k个与待处理像素的灰度差最小的像素，将这k个像素的灰度均值替代原来的像素值，譬如下图，给定3*3模板，k=5 需要注意的是，在模板上寻找与待处理像素最接近的k个像素点的时候，待处理像素也被包含在内，这是正确的，千万不要将其排除在外

In [15]:
t=conv(np.rollaxis(lena_noise,2),np.ones((3,3)),ratio=None) #平滑均值滤波
t=conv(t,np.ones((3,3)),ratio=None)
t=conv(t,np.ones((3,3)),ratio=None)
t=conv(t,np.ones((3,3)),ratio=None)
t=conv(t,np.ones((3,3)),ratio=None)
t=conv(t,np.ones((3,3)),ratio=None)
plt.figure(figsize=(15,15))
plt.subplot(121)
plt.imshow(t,cmap='gray')
t=cv2.cvtColor(lena_noise,cv2.COLOR_BGR2GRAY)
t=knn_smooth(t) #K近邻平滑滤波
t=knn_smooth(t)
t=knn_smooth(t)
t=knn_smooth(t)
t=knn_smooth(t)
t=knn_smooth(t)
plt.subplot(122)
plt.imshow(t,cmap='gray')
Out[15]:
<Figure size 1080x1080 with 0 Axes>
Out[15]:
<matplotlib.axes._subplots.AxesSubplot at 0x16773b204e0>
Out[15]:
<matplotlib.image.AxesImage at 0x16773b5eef0>
Out[15]:
<matplotlib.axes._subplots.AxesSubplot at 0x16773b5ef28>
Out[15]:
<matplotlib.image.AxesImage at 0x16773b642e8>

上述实验中，我对噪声图分别进行了6次均值滤波和6次K近邻均值滤波，因为只执行一次看不出模糊和不模糊的区别，多执行几次，效果就明显了

对称近邻平滑滤波器

image.png

最小方差平滑滤波器

设有一个模板，如果模板中的像素属于同一个区域，则模板中不包含边界像素，可以进行平滑处理；如果模板中的像素属于至少两个不同的区域，则模板中包含有边界像素，这时要对其进行保持，不进行平滑处理。要判断模板中的像素是否属于同一区域，一个最常用的方法是计算模板中所有像素的灰度方差，如果方差大，则表明模板像素属于不同区域的可能性大。考虑到景物边界的不规则性，选9个不同形状的模板。对9个模板所覆盖区域中的像素，分别计算其灰度分布方差，然后选择出方差为最小的模板中的像素灰度均值代替原像素值。下图展示了9个不同形状的平滑处理模板： image.png 其中"O"包裹的像素点为当前待处理像素点^

In [4]:
img=cv2.cvtColor(lena_noise,cv2.COLOR_BGR2GRAY)
t=knn_smooth(img) #K近邻平滑滤波
plt.figure(figsize=(15,15))
plt.subplot(121)
plt.imshow(t,cmap='gray')
t=min_gray_var(img,accelerate=True) #灰度最小方差平滑滤波
plt.subplot(122)
plt.imshow(t,cmap='gray')
Out[4]:
<Figure size 1080x1080 with 0 Axes>
Out[4]:
<matplotlib.axes._subplots.AxesSubplot at 0x2d7414cc550>
Out[4]:
<matplotlib.image.AxesImage at 0x2d7414b1c18>
Out[4]:
<matplotlib.axes._subplots.AxesSubplot at 0x2d741384080>
Out[4]:
<matplotlib.image.AxesImage at 0x2d741515940>

图像某一区域内灰度的方差确实能从一个角度反应区域内图像的对比度，方差越大，意味着图像区域内像素灰度变化越剧烈（灰度的平均值表示图像的总体亮度）

Sigma平滑滤波器

根据统计数学的原理，属于同一类别的元素的置信区间为均值附近$\pm2\sigma$范围。Sigma滤波器是构造一个模板，计算模板的标准差，置信区间为当前像素值的$\pm2\sigma$范围。将模板中落在置信区间范围内得像素的均值替换原来的像素值 image.png

In [8]:
img=cv2.cvtColor(lena_noise,cv2.COLOR_BGR2GRAY)
t=sigma_filter(sigma_filter(sigma_filter(img))) #执行了三次sigma平滑过滤，一次根本看不出有多大改变
plt.figure(figsize=(8,8))
plt.imshow(t,cmap='gray')
Out[8]:
<Figure size 576x576 with 0 Axes>
Out[8]:
<matplotlib.image.AxesImage at 0x194732b80f0>

小结：边界保持类平滑滤波器的核心是，尽可能地避免将平滑处理用于两个或多个不同区域，否则会导致边界模糊。可以采用不同形状结构判别，就我个人所做的实验来说，最小方差平滑滤波器和Sigma平滑滤波器的效果还是比较好的，但是计算复杂度也较高，所幸可以使用矩阵计算，但是sigma平滑滤波就难以使用矩阵计算

示例

In [21]:
demo=cv2.imread('./DB/image/22.jpg')
plt.figure(figsize=(15,15))
plt.subplot(121)
plt.imshow(demo[...,::-1])
plt.subplot(122)
ret=None
for i in range(3):
    t=sigma_filter(demo[:,:,i])[:,:,None]
    if i==0:
        ret=t
    else:
        ret=np.concatenate((ret,t),axis=2)
ret=media_filter(np.rollaxis(ret,2))
ret=np.rollaxis(ret,0,3)
plt.imshow(ret[...,::-1])
Out[21]:
<Figure size 1080x1080 with 0 Axes>
Out[21]:
<matplotlib.axes._subplots.AxesSubplot at 0x1f6a52434a8>
Out[21]:
<matplotlib.image.AxesImage at 0x1f6a530bac8>
Out[21]:
<matplotlib.axes._subplots.AxesSubplot at 0x1f6a47279b0>
Out[21]:
<matplotlib.image.AxesImage at 0x1f6a46b3e48>

尼玛这效果是怎么来的？ image.png

参考：
[1] https://blog.csdn.net/u012936765/article/details/53200918 生成噪声和各种滤波器
[2] http://www.ruanyifeng.com/blog/2017/12/image-and-wave-filters.html 图像与波形